source('00_util_scripts/mod_bplot.R')
source('00_util_scripts/mod_seurat.R')

proj.nm <- 'mission/ADRB2-FOH/'

sobj <- read_rds('mission/FPP/xiangya_sle_scRNA/data/human_SLE_derm_epi.rds')

sobj_epi <- read_rds('mission/FPP/xiangya_sle_scRNA/data/human_SLE_epi.rds')

sobj_epi |> DotPlot('ADRB2', cols = 'RdYlBu')

sobj_epi$sle |> unique()

sobj_epi |> VlnPlot('ADRB2', group.by = 'manual_main', pt.size = 0)

# All Adr genes -------
adr.hs.gene <- sobj_epi |> rownames() |> str_subset('^ADR[A-C]') |>
  str_sort()

sobj |> filter(manual_main != 'B_cell') |>
  mutate(manual_main = fct_relevel(manual_main, 'Keratinocytes')) |>
  DotPlot(adr.hs.gene, cols = 'RdYlBu', group.by = 'manual_main') +
  labs(x = 'Gene', y = 'Cell type',
       title = 'Adrenoceptors in human skin')

adr.hskin.expr <- last_plot() |>
  pluck('data') |>
  as_tibble()

adr.hskin.expr |>
  BubblePlot() +
  labs(title = 'Adrenoceptors in human skin') +
  theme_bw(base_size = 6) +
  RotatedAxis() +
  theme(legend.key.size = unit(2, 'mm'), legend.position = 'bottom')

proj.nm <- 'mission/ADRB2-FOH/'

publish_source_plot('human.skin.adr.genes.dotplot')

sobj_epi |>
  filter(manual_main != 'B_cell') |>
  mutate(manual_main = fct_relevel(manual_main, 'Keratinocytes')) |>
  bill.violin('ADRB2', group.by = manual_main) +
  labs(x = 'Cell type', y = 'Normalized expression',
       title = 'Human epidermis ADRB2 expression') +
  NoLegend()

## Dotplot in KC -------
sobj_keras |>
  mutate(sle == 'HC') |>
  DotPlot(c('ADRB2', late.kc, 'KRT5', 'KRT14'),
          cols = 'RdYlBu', cluster.idents = T) +
  labs(x = 'Gene', y = 'Keratinocyte clusters',
       title = 'ADRB2 and late KC markers in human skin KC')

kc.adrb2.diffmark <- last_plot() |>
  pluck('data')

# TODO: rewrite in tidyplots
kc.adrb2.diffmark |>
  filter(features.plot == 'KRT1') |>
  mutate(order = fct_reorder(id, avg.exp), id, .keep = 'none') |>
  right_join(kc.adrb2.diffmark) |>
  mutate(id = order) |>
  BubblePlot() +
  labs(title = 'ADRB2 & KC differentiation marker in human KC') +
  theme_jpub +
  RotatedAxis()

publish_source_plot('human.kc.adrb2.kcdiff.dotplot')

## define ADRB2-hi KC ----------
adrb2.exp <-
  read_csv('mission/ADRB2-FOH/results/human.kc.adrb2.kcdiff.dotplot.csv')

adrb2.exp |>
  filter(features.plot == 'ADRB2') |>
  ggplot(aes(avg.exp, avg.exp)) +
  geom_point() +
  geom_text(aes(label = id), nudge_y = .1)

## DEGA ----------
b2h.kc.deg <- sobj_keras |>
  FindMarkersAcrossVar(split.by = 'sle', group.by = 'seurat_clusters',
                       ident.1 = c(4,7,5))

b2h.kc.deg |>
  write_csv('mission/ADRB2-FOH/results/human.b2hvl.kc.deg.csv')

b2h.kc.deg |>
  filter(cluster == 'HC') |>
  mutate(p_val_adj = ifelse(p_val_adj == 0, 1e-300, p_val_adj)) |>
  plot_bill_volc(group1 = 'ADRB2-high KC',
                 group2 = 'ADRB2-low KC')

library(clusterProfiler)
hc.b2h.gsego <- b2h.kc.deg |>
  filter(cluster == 'HC', p_val_adj < .05) |>
  pull(avg_log2FC, name = gene) |>
  sort(decreasing = T) |>
  gseGO(ont = 'ALL', OrgDb = 'org.Hs.eg.db', keyType = 'SYMBOL')

hc.b2h.gsego <- hc.b2h.gsego |>
  simplify()

hc.b2h.gsego@result |>
  write_csv('mission/ADRB2-FOH/results/human.b2h.kc.gsego.csv')

hc.b2h.gsego@result |>
  plot_enrichment(metric = NES, base_col = 'blue') +
  labs(x = 'NES',
       title = 'Downregulated GO pathways\nADRB2-high KC vs ADRB2-low KC') +
  theme_bw()

hc.b2h.gsego@result |>
  filter(NES > 0, ONTOLOGY == 'BP') |>
  plot_enrichment(metric = NES) +
  labs(x = 'NES',
       title = 'Upregulated GO pathways\nADRB2-high KC vs ADRB2-low KC')

# GPCR in skin KC -----------
## human ----------
sobj <- read_rds('mission/FPP/xiangya_sle_scRNA/data/human_SLE_derm_epi.rds')

hpca <- celldex::HumanPrimaryCellAtlasData()

gpcrs <- query_uniprot_keyword('KW-0297')

sobj$sle |> unique()

hs_kc_deg <- sobj |>
  filter(sle == 'HC') |>
  FindMarkers(group.by = 'manual_main', ident.1 = 'Keratinocytes') |>
  as_tibble(rownames = 'gene')

hs_kc_mark_gpcr <- hs_kc_deg |>
  filter(gene %in% gpcrs$symbol, p_val_adj < .05, avg_log2FC > 1)

sobj |>
  filter(sle == 'HC') |>
  DotPlot(head(hs_kc_mark_gpcr$gene, n = 20), 
          cols = 'RdYlBu', group.by = 'manual_main') +
  RotatedAxis()

last_plot() |>
  pluck('data') |>
  mutate(id = fct_relevel(id, 'Keratinocytes', after = Inf)) |>
  BubblePlot(size = 7) +
  RotatedAxis() +
  labs(title = 'Top 20 GPCR specificly expressed in HC human skin')

hs_bc_deg <- sobj |>
  filter(sle == 'HC') |>
  FindMarkers(group.by = 'manual_main', ident.1 = 'B_cell') |>
  as_tibble(rownames = 'gene')

hs_bc_mark_gpcr <- hs_bc_deg |>
  filter(gene %in% gpcrs$symbol, p_val_adj < .05, avg_log2FC > 1)

sobj |>
  filter(sle == 'HC') |>
  DotPlot(head(hs_bc_mark_gpcr$gene, n = 10), 
          cols = 'RdYlBu', group.by = 'manual_main') +
  RotatedAxis()

last_plot() |>
  pluck('data') |>
  mutate(id = fct_relevel(id, 'B_cell', after = Inf)) |>
  BubblePlot(size = 7) +
  RotatedAxis() +
  labs(title = 'Top 10 GPCR specificly expressed in HC human skin B cells')

histamine_recp <- c('HRH1','HRH2','HRH3','HRH4')

sobj |>
  filter(sle == 'HC') |>
  DotPlot(histamine_recp, cols = 'RdYlBu', group.by = 'manual_main') +
  RotatedAxis()

last_plot() |>
  pluck('data') |>
  BubblePlot(size = 7) +
  RotatedAxis() +
  labs(title = 'Histamine receptors in HC human skin')

### HPCR ------------
hpca_gpcr <- hpca |>
  filter(.feature %in% gpcrs$symbol)

mean_gpcr <- hpca_gpcr |>
  summarise(mean_logcount = logtpm.mean(logcounts), .by = c(.feature, label.main))

bc_gpcr <- mean_gpcr |>
  slice_max(mean_logcount, by = .feature) |>
  filter(label.main == 'B_cell')

zs_gpcr <- mean_gpcr |>
  filter(.feature %in% bc_gpcr$.feature) |>
  mutate(z_score = scale(mean_logcount)[,1], .by = .feature)

zs_gpcr |>
  mutate(label.main = fct_relevel(label.main, 'B_cell')) |>
  ggplot(aes(label.main, .feature, fill = z_score)) +
  geom_tile() +
  scale_fill_distiller(palette = 'RdYlBu', values = pretty_distiller(zs_gpcr$z_score)) +
  rotate_x_text() +
  labs(x = 'Cell type', y = 'Gene', fill = 'Mean expr',
       title = 'Top expressed GPCRs in human B cells',
       subtitle = 'Human protein atlas')

## mouse -----------
gpcrs_mm <- query_uniprot_keyword('KW-0297', species = 'mouse')

sobj <- read_rds('mission/FPP/zww_sa_mice/data/aureus_mice_skin.rds')

sobj$orig.ident |> unique()

pbs_kc_deg <- sobj |>
  filter(orig.ident == 'PBS') |>
  FindMarkers(group.by = 'manual_main', ident.1 = 'Keratinocytes') |>
  as_tibble(rownames = 'gene')

pbs_kc_mrk_gpcr <- pbs_kc_deg |>
  filter(gene %in% gpcrs_mm$symbol, p_val_adj < .05, avg_log2FC > 1)

pbs_kc_mrk_gpcr |>
  head(n = 20) |>
  arrange(pct.2)
  filter(gene == 'Adgrl2')

sobj |>
  filter(orig.ident == 'PBS') |>
  DotPlot(head(pbs_kc_mrk_gpcr$gene, n = 20), 
          cols = 'RdYlBu', group.by = 'manual_main') +
  RotatedAxis()

last_plot() |>
  pluck('data') |>
  mutate(id = fct_relevel(id, 'Keratinocytes', after = Inf)) |>
  BubblePlot(size = 7) +
  RotatedAxis() +
  labs(title = 'Top 20 GPCR specificly expressed in PBS mouse skin')

sobj |>
  filter(orig.ident == 'PBS') |>
  DotPlot(str_to_title(histamine_recp),
          cols = 'RdYlBu', group.by = 'manual_main') +
  RotatedAxis()

last_plot() |>
  pluck('data') |>
  BubblePlot(size = 7) +
  RotatedAxis() +
  labs(title = 'Histamine receptors in PBS mouse skin')
